# TABLE 2
source('scripts/createDataset.R')

# ALIGNMENT
fileName <- "results/MANUSCRIPT/table_2_1.html"
add_vars <- c("partySPD", "partyGRUNE", "partyFDP", "partyLINKE", "partyAFD")
selectID <- df %>% 
  dplyr::mutate(ID = 1:nrow(.)) %>%
  dplyr::select(all_of(c(dep, run, contVarBase, add_vars)), ID) %>%
  subset(., complete.cases(.)) %>%
  pull(ID)
temp <- df %>% 
  mutate(ID = 1:nrow(.)) %>%
  subset(., ID %in% selectID) %>%
  select(all_of(c(dep, run, contVarBase, 'income2', add_vars))) %>%
  as.data.frame()
stargazer(temp, type = "html", summary.stat = c("n","min", 'max', 'mean', 'sd'), summary = T, out = fileName)

# TURNOUT
fileName <- "results/MANUSCRIPT/table_2_2.html"
temp <- df_turnout %>% select(turnout13) %>% as.data.frame()
stargazer(temp, type = "html", summary.stat = c("n","min", 'max', 'mean', 'sd'), summary = T, out = fileName)

############

# FIGURE 1
source('scripts/createDataset.R')

temp <- dat_temp %>%
  subset(., abs(x) <= bw) %>%
  mutate(x = floor(x/6)*6) %>%
  group_by(x) %>% summarize(y = mean(y, na.rm = T), n = n()) %>% ungroup() %>%
  mutate(D = case_when(x < 0 ~ 'not eligible for Mütterrente', T ~ 'eligible for Mütterrente'))
ggplot() +
  geom_vline(xintercept = 0) +
  geom_smooth(aes(x = (-1)*(x+1), y = y, color = D), size = 1, method = "lm",
              formula = y ~ poly(x, 2), se = F, level = 0.95, data = temp, alpha = 0.5) +
  geom_smooth(aes(x = (-1)*(x+1), y = y, color = D), size = 1, method = "lm", linetype = 'dashed',
              formula = y ~ poly(x, 1), se = F, level = 0.95, data = temp, alpha = 0.5) +
  geom_point(aes(x = (-1)*(x+1), y = y, color = D), size = 3, shape = 16, stroke = 0, data = temp) +
  scale_color_manual(values = c('eligible for Mütterrente' = 'black', 
                                'not eligible for Mütterrente' = 'gray50')) +
  scale_size_continuous(guide = "none") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 5L), breaks = seq(0, 1, 0.1), limits = c(0, 0.8)) +
  scale_x_continuous(breaks = seq(-480, 480, 12)) +
  xlab('number of months prior/after January 1992') +
  ylab('alignment with the \npledge-making party') +
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom')
ggsave("results/MANUSCRIPT/figure_1.png", width = 10, height = 4)

############

# FIGURE 2
source('scripts/createDataset.R')

years <- seq(1999, 2017, 1)
res <- data.frame(stringsAsFactors = F)
for(i in 1:length(years)){
  start <- years[i]
  end <- years[i]+2
  paste0('YEAR OF SURVEY: ', start) %>% print()
  df_temp <- df_dur %>% 
    subset(., abs(x) <= bw_cct) %>%
    subset(., syear %in% c(start:end)) %>%
    select(x, y, pid, hid, all_of(contVarBase)) %>%
    subset(., complete.cases(.)) %>%
    mutate(D = case_when(x < 0 ~ 'not eligible for Mütterrente (control group)', 
                         T ~ 'eligible for Mütterrente (treatment group)')) %>%
    group_by(D) %>% summarise_at("y",funs(mean, sd), na.rm = T) %>% ungroup() %>%
    mutate(PERIOD = paste(start, end, sep = '-'))
  res <- res %>% bind_rows(., df_temp)
}

# CREATE PLOT
res <- res %>% mutate(D = factor(D), PERIOD = factor(PERIOD))
ggplot() +
  geom_rect(aes(xmin = "2009-2011", xmax = "2012-2014", ymin = -Inf, ymax = Inf), fill = 'gray70', alpha = 0.4) +
  geom_line(data = res, aes(x = PERIOD, y = mean, color = D, group = D), size = 2) +
  geom_point(data = res, aes(x = PERIOD, y = mean, color = D, group = D), size = 3.5) +
  scale_color_manual(values = c('eligible for Mütterrente (treatment group)' = 'black', 
                                'not eligible for Mütterrente (control group)' = 'gray50'),
                     drop = F) +
  scale_y_continuous(breaks = seq(0, 1, 0.05), limits = c(0, 0.45), labels = scales::percent) +
  xlab('years of survey') + 
  ylab('alignment with the \npledge-making party') + 
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom',
                     axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("results/MANUSCRIPT/figure_2.png", width = 6, height = 4)

############

# FIGURE 3
source('scripts/createDataset.R')

# ALIGNMENT
reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg[[1]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 1, 
                                 covariates = contVarBase %>% paste(., collapse = '+'))
reg[[1]]$param <- 'baseline'
reg[[2]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 2, 
                                 covariates = contVarBase %>% paste(., collapse = '+'))
reg[[2]]$param <- 'quadratic form'
reg[[3]] <- dat_temp1 %>% rdd_reg_lm(slope = "same", bw = bw, order = 1, 
                                 covariates = contVarBase %>% paste(., collapse = '+'))
reg[[3]]$param <- 'equal slopes'
reg[[4]] <- dat_temp1 %>% rdd_reg_np(slope = "same", bw = bw, 
                                     covariates = contVarBase %>% paste(., collapse = '+'))
reg[[4]] <- reg[[4]][['RDDslot']][['model']]
reg[[4]]$param <- 'non-parametric'
reg[[5]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 1)
reg[[5]]$param <- 'w/o controls'
reg_lin_alignment <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  mutate(dep = 'alignment', model = 'linear probability model')

# TURNOUT
reg <- list()
df_turnout <- df_turnout %>% mutate(y = get('turnout13'))
dat_temp1 <- rdd_data(df_turnout$y, df_turnout$x, covar = df_turnout[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg[[1]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 1, 
                                     covariates = contVarBase %>% paste(., collapse = '+'))
reg[[1]]$param <- 'baseline'
reg[[2]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 2, 
                                     covariates = contVarBase %>% paste(., collapse = '+'))
reg[[2]]$param <- 'quadratic form'
reg[[3]] <- dat_temp1 %>% rdd_reg_lm(slope = "same", bw = bw, order = 1, 
                                     covariates = contVarBase %>% paste(., collapse = '+'))
reg[[3]]$param <- 'equal slopes'
reg[[4]] <- dat_temp1 %>% rdd_reg_np(slope = "same", bw = bw,
                                     covariates = contVarBase %>% paste(., collapse = '+'))
reg[[4]] <- reg[[4]][['RDDslot']][['model']]
reg[[4]]$param <- 'non-parametric'
reg[[5]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 1)
reg[[5]]$param <- 'w/o controls'
reg_lin_turnout <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  mutate(dep = ' turnout', model = 'linear probability model')

# CREATE PLOT
plot_df <- bind_rows(reg_lin_alignment, reg_lin_turnout) %>%
  mutate(param = factor(param, levels = c("non-parametric", "equal slopes", "quadratic form", "w/o controls", "baseline")))
p1 <- ggplot() +
  geom_vline(xintercept = c(0)) +
  geom_linerange(data = plot_df, aes(xmin = cil95, xmax = ciu95, y = dep, group = param), 
                 size = 1.5, position = position_dodge(width = 0.5)) +
  geom_linerange(data = plot_df, aes(xmin = cil99, xmax = ciu99, y = dep, group = param), 
                 size = 0.75, position = position_dodge(width = 0.5)) +
  geom_point(data = plot_df, aes(x = estimate, y = dep, shape = param), 
             size = 3, stroke = 1.5, position = position_dodge(width = 0.5), fill = 'white') +
  scale_shape_manual(values = c(21:24, 15), drop = F) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1), breaks = seq(-0.5, 0.5, 0.05), limits = c(-0.15, 0.25)) +
  xlab('treatment effect') +
  ylab('') +
  guides(col = guide_legend(nrow = 1, byrow = TRUE, reverse = T),
         shape = guide_legend(nrow = 1, byrow = TRUE, reverse = T)) +
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom')
fileName_plot <- paste0('results/MANUSCRIPT/figure_3.png')
ggsave(fileName_plot, p1, width = 7, height = 3.5)

############

# FIGURE 4
source('scripts/createDataset.R')

df1 <- df
res <- data.frame(stringsAsFactors = F)
parties <- c("partyUNION", "partyGRUNE", "partySPD", "partyFDP", "partyLINKE", 'partyAFD')
for(i in 1:length(parties)){
  
  df <- df1 %>% mutate(y = get(parties[i]))
  reg <- list()
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  
  ###
  
  reg[[1]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 1,
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[1]]$param <- '1st | separate slopes | parametric (baseline)'
  reg[[2]] <- dat_temp1 %>% rdd_reg_lm(slope = "separate", bw = bw, order = 2,
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[2]]$param <- '2nd | separate slopes | parametric'
  reg[[3]] <- dat_temp1 %>% rdd_reg_lm(slope = "same", bw = bw, order = 1,
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[3]]$param <- '1st | same slopes | parametric'
  reg[[4]] <- dat_temp1 %>% rdd_reg_np(slope = "separate", bw = bw,
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[4]] <- reg[[4]][['RDDslot']][['model']]
  reg[[4]]$param <- '1st | separate slopes | nonparametric'
  
  temp <- reg %>% 
    extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
    mutate(PARTY = gsub('party', '', parties[i]), model = 'linear probability model')
  res <- bind_rows(res, temp) %>% mutate(DEP = 'alignment')
}

temp <- res %>% 
  subset(., !grepl('2nd', param)) %>% 
  subset(., grepl('baseline', param)) %>%
  mutate(PARTY = PARTY %>% factor(., levels = c(gsub('party', '', parties)))) %>%
  mutate(PARTYgroup = factor(PARTY, levels = c("UNION", "GRUNE", "SPD", "FDP", "LINKE", "AFD")))
p1 <- ggplot() +
  geom_vline(xintercept = 0) +
  geom_errorbar(data = temp, aes(y = reorder(PARTY, desc(PARTY)), x = estimate, 
                                 xmin = cil95, xmax = ciu95, color = PARTYgroup), size = 1.5, width = 0) +
  geom_errorbar(data = temp, aes(y = reorder(PARTY, desc(PARTY)), x = estimate, 
                                 xmin = cil99, xmax = ciu99, color = PARTYgroup), size = 0.75, width = 0) +
  geom_point(data = temp, aes(y = reorder(PARTY, desc(PARTY)), x = estimate, shape = PARTYgroup, color = PARTYgroup), fill = 'white', size = 4, stroke = 1.5) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1), breaks = seq(-0.5, 0.5, 0.05), limits = c(-0.25, 0.25)) +
  scale_y_discrete(labels = c('UNION' = 'conservatives\n[pledge-maker]', 'GRUNE' = 'greens', 'SPD' = 'social\ndemocrats', 
                              'FDP' = 'liberals', 'LINKE' = 'left socialists', 'AFD' = 'nationalists'), 
                   limits = c('AFD', 'LINKE', 'FDP', 'SPD', 'GRUNE', 'UNION'),
                   drop = FALSE) +
  xlab('treatment effect') +
  ylab('party') +
  scale_color_manual(drop = FALSE,
                     values = c('UNION' = '#000000', 'GRUNE' = '#46962b', 'SPD' = '#EB001F',
                                'FDP' = '#F6BE00',  'LINKE' = '#BE3075', 'AFD' = '#009EE0'),
                     labels = c('UNION' = 'CDU/CSU (conservatives) [pledge-maker]', 'GRUNE' = 'Die Grünen (greens)', 'SPD' = 'SPD (social democrats)',
                                'FDP' = 'FDP (liberals)', 'LINKE' = 'Die Linke (left socialists)', 'AFD' = 'AfD (nationalists/far-right)')) +
  scale_shape_manual(drop = FALSE,
                     values = c(15, 21:25),
                     labels = c('UNION' = 'CDU/CSU (conservatives) [pledge-maker]', 
                                'GRUNE' = 'Die Grünen (greens)', 'SPD' = 'SPD (social democrats)', 
                                'FDP' = 'FDP (liberals)', 'LINKE' = 'Die Linke (left socialists)', 
                                'AFD' = 'AfD (nationalists/far-right)')) +
  guides(shape = guide_legend(nrow = 2, byrow = TRUE, drop = F)) +
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom')
fileName_plot <- paste0('results/MANUSCRIPT/figure_4.png')
ggsave(fileName_plot, p1, width = 7.5, height = 3.5)

############

# FIGURE 5
source('scripts/createDataset.R')

mediatorvec <- c("age_olderthan47", "nchild_lessthan2", "nchild_morethan1", 
                 "married_married", "married_notmarried", "income2_lessthan850")
contVar <- contVarBase
robust_temp <- 'HC1'
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp) %>% 
    summary() %>% .[["coefficients"]] %>% as.data.frame() %>% 
    purrr::set_names('estimate', 'se', 'tval', 'pval') %>% 
    subset(., row.names(.) == 'D') %>%
    mutate(mediator = mediator_temp)
}
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg[[length(reg)+1]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp) %>% 
  summary() %>% .[["coefficients"]] %>% as.data.frame() %>% 
  purrr::set_names('estimate', 'se', 'tval', 'pval') %>% 
  subset(., row.names(.) == 'D') %>%
  mutate(mediator = 'baseline')
plot_df <- bind_rows(reg) %>%
  mutate(cil95 = estimate-(1.96*se), ciu95 = estimate+(1.96*se),
         cil99 = estimate-(2.576*se), ciu99 = estimate+(2.576*se)) %>%
  mutate(param = mediator) %>%
  mutate(vartype = case_when(param == "nchild_lessthan2" ~ 'family size', 
                             param == "nchild_morethan1" ~ 'family size',
                             param == "married_married" ~ 'marital status',
                             param == "married_notmarried" ~ 'marital status', 
                             param == "income2_lessthan850" ~ 'income level', 
                             param == "baseline" ~ 'baseline')) %>%
  mutate(param = case_when(param == "nchild_lessthan2" ~ 'single-child family', 
                           param == "nchild_morethan1" ~ 'multi-child family',
                           param == "married_married" ~ 'married',
                           param == "married_notmarried" ~ 'single/divorced/widowed', 
                           param == "income2_lessthan850" ~ 'less than 850 Euro/month', 
                           param == "baseline" ~ 'baseline')) %>%
  mutate(param = factor(param, levels = c("baseline", 
                                          "less than 850 Euro/month", 
                                          "single/divorced/widowed", 
                                          "multi-child family", 
                                          "married",
                                          "single-child family"))) %>%
  subset(., vartype != 'age')

vartypes <- list(
  c("baseline"),
  c("income level", "baseline"),
  c("marital status", "income level", "baseline"),
  c("family size", "marital status", "income level", "baseline")
)
temp <- plot_df %>%
  mutate(vartype = factor(vartype, levels = c("family size", "marital status", "income level", "baseline")))
p1 <- ggplot() +
  geom_vline(xintercept = c(0)) +
  geom_linerange(data = temp, aes(xmin = cil95, xmax = ciu95, y = param, group = param, color = vartype), 
                 size = 1.5, position = position_dodge(width = 0.5)) +
  geom_linerange(data = temp, aes(xmin = cil99, xmax = ciu99, y = param, group = param, color = vartype), 
                 size = 0.75, position = position_dodge(width = 0.5)) +
  geom_point(data = temp, aes(x = estimate, y = param, color = vartype), 
             size = 3, stroke = 1.5, position = position_dodge(width = 0.5), fill = 'white') +
  scale_color_jama(drop = F) +
  scale_shape_manual(values = c(21:24, 15)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1), breaks = seq(-0.5, 0.5, 0.05),
                     limits = c(-0.15, 0.5)) +
  xlab('treatment effect') +
  ylab('sub-sample') +
  scale_y_discrete(limits = rev, drop = F) +
  guides(col = guide_legend(nrow = 1, byrow = TRUE, reverse = T),
         shape = guide_legend(nrow = 1, byrow = TRUE, reverse = T)) +
  theme_bw() + theme(legend.title = element_blank(), legend.position = 'bottom')
fileName_plot <- paste('results/MANUSCRIPT/figure_5.png')
ggsave(fileName_plot, p1, width = 7, height = 3.5)

############

# FIGURE 6
source('scripts/createDataset.R')

# DEFINE START AND END DATE
reg <- list()
for(j in 1:length(years)){
  dat_temp1 <- defineDatasetDuration()
  reg[[length(reg)+1]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[length(reg)]]$param <- paste(years[j], years[j]+2, sep = "-")
}
res_union <- extractRegResultsForPlot(reg, robust_temp = robust_temp, exponentiate = F) %>%
  mutate(param = param %>% factor(param), PARTY = 'UNION')

##

reg <- list()
for(j in 1:length(years)){
  df_dur <- df_dur %>% mutate(y = partyGRUNE)
  dat_temp1 <- defineDatasetDuration()
  reg[[length(reg)+1]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[length(reg)]]$param <- paste(years[j], years[j]+2, sep = "-")
}
res_grune <- extractRegResultsForPlot(reg, robust_temp = robust_temp, exponentiate = F) %>%
  mutate(param = param %>% factor(param), PARTY = 'GRUNE')

temp <- bind_rows(res_union, res_grune)

p <- ggplot() +
  geom_rect(aes(xmin = "2009-2011", xmax = "2012-2014", ymin = -Inf, ymax = Inf), fill = 'gray70', alpha = 0.4) +
  geom_hline(yintercept = 0) +
  geom_errorbar(data = temp, aes(x = param, y = estimate, ymin = cil95, ymax = ciu95, 
                                 group = PARTY, color = PARTY), size = 1.5, width = 0) +
  geom_errorbar(data = temp, aes(x = param, y = estimate, ymin = cil99, ymax = ciu99,
                                 group = PARTY, color = PARTY), size = 0.75, width = 0) +
  geom_point(data = temp, aes(x = param, y = estimate, group = PARTY, color = PARTY), size = 3) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), 
                     breaks = seq(-1, 1, 0.05), limits = c(-0.2, 0.2)) +
  scale_x_discrete(drop = F) +
  scale_color_manual(values = c('GRUNE' = '#46962b', 'UNION' = '#000000'), 
                     labels = c('UNION' = 'CDU/CSU (conservatives) [pledge-maker]',
                                'GRUNE' = 'Die Grünen (greens)'),
                     drop = F) +
  xlab('survey period') + ylab('treatment effect') + 
  theme_bw() +
  theme(legend.position = 'bottom', legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
fileName_plot <- paste0('results/MANUSCRIPT/figure_6.png')
ggsave(fileName_plot, plot = p, width = 7, height = 4)

